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A generalization of the 3 D Euler-Voigt-a model is obtained by introducing derivatives of arbitrary 
order /? (instead of 2) in the Helmholtz operator. The /? —> 00 limit is shown to correspond to 
Galerkin truncation of the Euler equation. Direct numerical simulations (DNS) of the model are 
performed with resolutions up to 2048 3 and Taylor-Green initial data. DNS performed at large /3 
demonstrate that this simple classical hydrodynamical model presents a self-truncation behavior, 
similar to that previously observed for the Gross-Pitaeveskii equation in Krstulovic and Brachet 
[Phys. Rev. Lett. 106, 115303 (2011)]. The self-truncation regime of the generalized model is 
shown to reproduce the behavior of the truncated Euler equation demonstrated in Cichowlas et al. 
[Phys. Rev. Lett. 95, 264502 (2005)]. The long-time growth of the self-truncation wavenumber k s t 
appears to be self-similar. 

Two related a-Voigt versions of the EDQNM model and the Leith model are introduced. These 
simplified theoretical models are shown to reasonably reproduce intermediate time DNS results. 
The values of the self-similar exponents of these models are found analytically. 


I. INTRODUCTION 

Classical Galerkin-truncated systems have been stud¬ 
ied since the early 50’s in fluid mechanics. In this context, 
the (time reversible) Euler equation describing spatially- 
periodic classical ideal fluids is known to admit, when 
spectrally truncated at wavenumber £; max , absolute equi¬ 
librium solutions with Gaussian statistics and equipar- 
tition of kinetic energy among all Fourier modes [1-4]. 
Furthermore, the dynamics of convergence toward equi¬ 
librium involves a direct energy cascade toward small- 
scales and contains (long-lasting) transient that mimic 
(irreversible) viscous effects that are produced by the 
“gas” of high-wavenumber partially-thermalized Fourier 
modes generating (pseudo) dissipative effects [5-7]. 

In the case of superfluids, the relevant equation is 
the so-called truncated (or Galerkin-projected) Gross- 
Pitaevskii equation (TGPE). In the TGPE case, abso¬ 
lute equilibrium can also be obtained by a direct energy 
cascade, in a way similar to that of the truncated Euler 
case, with final thermalization accompanied by vortex 
annihilation. Furthermore, increasing the amount of dis¬ 
persion produces a slowdown of the energy transfer at 
small-scales inducing a bottleneck and a partial thermal¬ 
ization that is independent of the truncation wavenum¬ 
ber and takes place below a ‘self-truncation’ wavenumber 
k s t{t ) that is observed to slowly increase with time [8-10]. 

The purpose of the present paper is to find and study 
such self-truncation phenomena in the simpler context 


of classical hydrodynamics of an ideal fluid. This is ob¬ 
tained by using equation of motion of the Euler type. To 
wit, we study here a simple generalization of the standard 
3D Euler-Voigt-a (non dissipative) model [11, 12]. 

Note that various dissipative Navier-Stokes-Voigt-a 
regularizations have been proposed in the last decade as 
efficient subgrid-scale models in order to address the clas¬ 
sical turbulence closure problem, both in hydrodynamics 
[13] and in magnetohydrodynamics [14]. Results on the 
regularity properties of this type of dissipative models 
are given in [15] and the statistical solutions are shown 
to converge to those of Navier-Stokes (when a —> 0) in 
[16, 17]. For large values of a, dissipative Voigt models 
are known to inhibit and reduce the transfer of energy to 
the small scales. 

Compared to the Euler equations, the conservative 3D 
Euler-Voigt-a model also penalizes the formation of small 
scales. We show that this penalization is enough to pro¬ 
duce a self truncation regime. Our main findings are 
that the self-truncation regime of this generalized model 
reproduces the behavior of the truncated Euler equation 
[5]. The long-time behavior of the energy spectrum ap¬ 
pears to be self-similar. 

To understand this self-similarity we further intro¬ 
duce two different models that are a-Voigt versions of 
the Eddy-Damped Quasi-Normal Markovian (EDQNM) 
model [4, 18] and the Leith model [19], respectively. Both 
models are shown to present behaviors that are similar 
to that of the Euler-Voigt-a model. The relative simplic- 


2 


ity of these models allows us to determine the analytical 
values of the self-similar exponents. 

The paper is organized as follows. Section II is de¬ 
voted to our generalized model. Basic definitions of the 
Euler-Voigt-a model are given in section II A. Numeri¬ 
cal methods and performed computations are detailed in 
section IIB. Our results on the self truncation regime 
are described in section IIC. The long-time behavior is 
studied in section IID. Related theoretical models are 
presented in section III. Section III A is devoted to the 
aV-EDQNM model and section IIIB to the aV-Leith 
model. Finally, our main results are summarized in sec¬ 
tion IV where we give our conclusions. 


II. EULER-VOIGT-a MODEL 


A. Definition of the model 


The standard 3 D Euler-Voigt-a model [11, 12] is a 
partial differential equation for the 3 D velocity field 
u (x,y,z,t) that explicitly reads 

r) 11 

(l-a 2 V 2 )^ = -(u-V)u-Vp 

V • u = 0. (1) 

The operator — a 2 V 2 (as we will see later) suppresses 
the formation of scales smaller than a. The associated 
wave number to this scales is denoted k a = a -1 . We 
refer to the operator in Eq.l as the a-term. 

Let us now define the generalized 3 D Euler-Voigt-a 
model: 

(1 + (-a 2 V 2 )^ = -(u • V)u - Vp 

V • u = 0, (2) 


where the power /? is an even integer. We refer to /3 as 
the penalization exponent as its increase enhances the 
suppression of small scale generation. When a = 0, the 
generalized model Eq. (2) reduces to the standard 3D 
incompressible Euler equations 

<9u 

——|-uVu = — Vp, V • u = 0. (3) 

at 

We consider here spatially-periodic solutions defined 
in the domain 12 = [0,27r] 3 . The kinetic energy spec¬ 
trum E(k,t) associated to (3) is defined as the sum over 
spherical shells 


E o%t) = \ Y l S (M)| 2 , ( 4 ) 

kezs 

k—l/2<|k|<k+l/2 


and the energy 
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|u(x, t)\ 2 d 3 x 


l j Y 

kez 3 


is independent of time when u satisfies the 3D Euler 
equations (3). The conserved energy associated to the 
generalized Euler-Voigt-a model (2) is straightforward 
to obtain and reads in physical space 

Ecc = ]3 [ u • [1 + (-a 2 V 2 )^]ud 3 i, 

and in spectral space 

^ = ^E[ 1 + ( a ^]i a ( k ’ < )i 2 - ( 5 ) 

kez 3 

Consequently, the generalized energy spectrum is defined 
as 

E a (k,t) = i Y [1 + (afc) /3 ]|u(k,i)| 2 . (6) 

kez 3 

k-l/2<|k|<k+l/2 

In the following, we refer to E a (k,t) as the energy spec¬ 
trum and E(k,t) as the kinetic energy spectrum. Equa¬ 
tions (2) also conserve the generalized helicity 

H a = ^Y^ 1 + i ak ) P ] S ( k ’ *) • £(-k, t). (7) 


In this work we only consider flows with H a = 0. 

Let us remark that the differential operator multi¬ 
plying the r.h.s. of our generalized 3D Euler-Voigt-a 
model (2) can be written in Fourier space as 1 + (ak)& = 
1 + (k/k a )P. The formal limit /3 —> oo of Eq. (2) thus 
corresponds to a standard spherical Galerkin truncation 
(u(k) = 0 for |k| > fcmax) of the Euler equation (3) at 
fcmax = k a . Note that a somewhat similar generalization 
(but involving the use of a high powers of the Laplacian 
in the dissipative term of forced hydrodynamical equa¬ 
tions) have been studied in reference [20]. 

It is well known that the truncated Euler equation ad¬ 
mits statistically stationary solutions given by the mi- 
crocanonical distribution determined by the invariants 
[1, 3]. These solutions are the so called absolute equilib¬ 
rium and lead to equipartition of energy among Fourier 
modes. The Euler-Voigt-a model considered as a trun¬ 
cated system, also admits absolute equilibrium solu¬ 
tions. When fully thermalized, Fourier modes can be 
described by the canonical Gibbs distribution u(k) ~ 
Z~ x exp {— /3E a [u]}, where Z is the partition function 
[21]. As the invariant energy is quadratic, the absolute 
equilibrium is Gaussian and the Fourier modes are inde¬ 
pendent. This leads to the spectra 


E a {k) 


3 E a k 2 


or E(k) 


3 E a k 2 

fcmax(! + ^kP )' 


( 8 ) 


The large scale behavior of the kinetic energy spectrum 
thus depends on the value of the penalization exponent 
P as E(k) ~ k 2 ~ lj . Therefore, in thermal equilibrium the 
small scales of u (i.e. k k a ) are penalized when /3 > 0. 

Taking into account that, in Fourier space, the dif¬ 
ferential operator in (2) can be defined for real values of 
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/3 > 0, the choice /? = 11/3 yields an absolute equilibrium 
E(k) ~ fc -5 / 3 and thus a fully thcrmalized field following 
Kolmogorov scaling. In the same vein, choosing in two- 
dimensions f) = 2/3 also yields Kolmogorov scaling, this 
time for the equipartition of enstrophy. This can repre¬ 
sent an interesting alternative to the fractal decimation 
method that was used in reference [22]. 


B. Numerical method 

The generalized 3D Euler-Voigt-a equations (2) are 
solved numerically using standard [23] pseudo-spectral 
methods with resolution N. Time marching is performed 
using a second-order Runge-Kutta scheme and the solu¬ 
tions are spherically dealiased by suppressing, at each 
time step, the modes for which the wave-vector exceeds 
two-thirds of the maximum wave-number N/2 (thus a 
2048 3 run is truncated at |k| > fc max = 682,see e.g. [24] 
for details on the numerical method). 

We consider here solutions of Eq. (2) that correspond 
to the so-called Taylor-Green (TG) [25] (27r-periodic) ini¬ 
tial data u(x, y, z, 0) = u TG (x, y, z), with 

u TG = (sin(x) cos(y) cos(z), — cos(x) sin(y) cos(z), 0). 

(9) 

The simulations reported in this paper were performed 
using a special purpose symmetric parallel code devel¬ 
oped from that described in [24, 26-28]. The code uses 
the symmetries of the Taylor-Green initial data to speed¬ 
up computations and optimize memory usage. The work¬ 
load for a timestep is (roughly) twice that of a general pe¬ 
riodic code running at a quarter of the resolution. Specif¬ 
ically, at a given computational cost, the ratio of the 
largest to the smallest scale available to a computation 
with enforced Taylor-Green symmetries is enhanced by a 
factor of 4 in linear resolution. This leads to a factor of 32 
savings in total computational time and memory usage. 
The code is based on FFTW and a hybrid MPI-OpenMP 
scheme derived from that described in [29]. At resolution 
2048 3 we used 512 MPI processes, each process spawning 
8 OpenMP threads. 

When compared with standard Euler equation (3) runs 
that were performed in reference [24], the only com¬ 
putational advantage of the the generalized 3D Euler- 
Voigt-a model (2) stems from the much weaker Courant- 
Friedrichs-Lewy (CFL) condition Uk At < C on the time 
step At which is conditioned by k = k a = a” 1 rather 
than k = k m ax , when k a « k m ax . We have performed 
a number of high-resolution runs that are summarized in 
Table I. The CFL condition A tk a ~ 0.1 was enough to 
insure both stability and energy conservation up to 0.3% 
in the worst case. 


C. Self truncation 

A first indication on the dynamics of the generalized 
Euler-Voigt-a model (2) with TaylorGreen initial data 


Run 

Res. 

P 

ka 

^max 

Run 

Res. 

P 

ka 

^max 

Euler 

1024 

- 

- 

30 

11 

2048 

2 

100 

15 

1 

1024 

2 

20 

50 

12 

2048 

2 

200 

15 

2 

1024 

2 

40 

50 

13 

2048 

4 

50 

15 

3 

1024 

2 

80 

50 

14 

2048 

4 

100 

15 

4 

1024 

4 

20 

50 

15 

2048 

4 

200 

15 

5 

1024 

4 

40 

50 

16 

2048 

6 

50 

15 

6 

1024 

4 

80 

50 

17 

2048 

6 

100 

15 

7 

1024 

6 

20 

50 

18 

2048 

6 

200 

15 

8 

1024 

6 

40 

50 

19 

512 

2 

4 

2300 

9 

1024 

6 

80 

50 

20 

512 

4 

4 

2300 

10 

2048 

2 

50 

15 

21 

512 

6 

4 

2300 


TABLE I. List of runs of the generalized 3 D Euler-Voigt-a 
model (2) with TaylorGreen initial data (9) and maximum 
integration time t max . 


(9) is given by the behavior of the Energy spectra for a 
run at resolution 1024 3 with /3 = 4 and k a = 80 (run 6 
Table I) that is displayed in Fig.l. 

10 ' 

10 ' 

10 ' 

' a 

LLI 

10 " 

10 ' 


FIG. 1. (Color online) Temporal evolution (indicated by ar¬ 
rows) of the energy spectrum E a (k) for /3 = 4 and k a = 80. 
Resolution 1024 3 (fc max = 342). The dashed lines respectively 
display the Kolmogorov AC 5 / 3 and the equipartition k 2 scal¬ 
ing. The self-truncation wavenumber is indicated by the small 
vertical arrow. 

Different regimes are clearly observed. First the energy 
is transferred toward small scales as in the standard Euler 
equation evolution (t < 5). Then the energy reaches the 
wavenumber k a = 80 and the a-term in (2) starts to sup¬ 
press the energy transfer for k > k a . As a consequence, 
the energy piles up around this wavenumber, similarly to 
the truncated Euler case with a cut-off ~ k a . A compat¬ 
ible Kolmogorov Ac -5 / 3 scaling is observed at large scales 
(t = 15), followed by a partially thermalized zone in k 2 
extending up to k s t . The wavenumber fc st (f) then slowly 
grows from its initial value k a until it eventually reaches 
the simulation cut-off fc max . For t —> oo the system fully 
thermalises independently of the parameters (data not 
shown) and the spectrum is then described by the abso¬ 
lute equilibrium (8). It is conspicuous that for k > fc st 
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the energy quickly decays and the partial thermalization 
regime is thus independent of the simulation cut-off fc max . 
We thus refer to k st as the self-truncation wavenumber. 
A further discussion and justification for its name will be 
given later (see below, paragraph following Eq.10). 

Figure 2 displays the spectra of simulations at reso¬ 
lution 2048 3 for different values of /3 and k a taken at 
time t = 10.2 [30]. The self-truncation is apparent for 



FIG. 2. (Color online) Energy spectrum E a (k ) versus k for 
/3 = 2 and k a = 50, 100 and 200 at t = 10.2. (b) Same con¬ 
ditions than in (a) but /3 = 4. c) Same conditions than in (a) 
but /? = 6. The dashed-black lines corresponds to Kolmogorov 
scaling E a (k) ~ fc~ 5//3 and energy equipartition E a (k) ~ k 2 . 
The thermalization wave number kth and selftruncation k B t 
are indicated by arrows in (c). 

the smaller k a at all values of f3 and, for larger /3 > 4, 
at all values of k a . The three ranges delimited by the 
thermalization wavenumber fc t h and the self-truncation 
wavenumber fc st are clearly visible (see Fig.2.c) at this 
high resolution. Formally, these zones correspond to the 
Kolomogorov regime (fc -C fc t h), the thermalization range 
( fcth *C k <C k s t) and the exponential energy decay 
(k fc s t). It is remarkable that, at this intermediate 


time (t = 10.2), the /3 = 6 case appears to be already 
behaving somewhat like what is expected of the (3 = oo 
limit. Indeed, Fig. 2(c) is reminiscent of the previously 
studied truncated-Euler case (see Fig.l of reference [5]). 
Of course, in the truncated Euler case no third decreas¬ 
ing zone exists as, when /? = oo, E a (k) = 0 for k > k a . 
In order to make the comparison with the (3 = oo limit 
more quantitative, following [5], we have computed the 
thermalized energy E t g(t) = J2k th (t)<k E{k,t) and the 

effective dissipation eft) = • The time evolutions 

of fcth, E t h and e are presented on Fig.3. The agreement 



t t 


FIG. 3. (Color online) Temporal evolution of (a) k t h, (b) 
.Eth, (d) e = g ) h , (e) k t h/kg with estimated based on the 
self-truncation wavenumber k s t, see Eq.(10). 

with the truncated Euler data appears to be good (com¬ 
pare with Figs 2 and 4 in reference [5]). These quantities 
do not appreciably depend on the value of {3. 

In references [5, 31], an effective generalized Navier- 
Stokes model for the dissipative dynamics of modes k 
close to kth(t) was suggested for the original Euler case 
with fixed truncation at k = fc m ax- In this case, the 
effective viscosity of the model was given by v e g = 
VAth/fcmax- If we assume that the generalized large- (3 
case behaves similarly to the Euler case truncated at fc st , 
we find that the dissipative wavenumber kg should be 
given by 

fc d ~ e^^s/E^/kst)- 3 /* (10) 

The consistency of this estimation of the effective dissi¬ 
pation with the results displayed in Fig.3a-c requires that 
kg ~ fcth- The ratio fcth/fcd is displayed on Fig.3.d. It is 
indeed of order unity and reasonably constant in time. 
Thus the larg e-/3 dynamics of Eq. (2) is seen to emulate 
the dynamics of the Euler equation, spectrally-truncated 
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FIG. 4. (Color online) (a) Exponential decay of E a (k) at 
t = 7.8. b) Time Evolution of energy spectrum fit parameter 5 
of Eq.(ll): Horizontal lines correspond to 5k max = 2 (dashed 
black line) and 5k a = 2 (dot-dashed black line). Exponential 
law (dot black line) 5 = 2.7 exp(—1/0.56) from reference [24]. 

at fcmax = k s t- For this reason, we call the regime where 
this behavior takes place the self-truncation regime. 

In the self-truncation regime the energy spectrum 
E a (k) sharply decreases for for k > k st (see Fig.2). In or¬ 
der to quantify this behavior we have used the so-called 
analyticity strip (AS) method [32]. Let us recall that 
the idea is to monitor the ‘width of the analyticity strip’ 
<5(> 0) as a function of time, effectively measuring a ‘dis¬ 
tance to the singularity’ [33]. Using spectral methods 

[23] , S(t) is obtained directly from the high-wavenumber 
exponential fall off of the spatial Fourier transform of the 
solution [34]. The AS method has been used to numeri¬ 
cally study putative Euler singularities (see e.g. reference 

[24] for implementation details). The common procedure 
is to perform a least-square fit at each time t on the log¬ 
arithm of the energy spectrum E a (k,t), using the func¬ 
tional form 

In E a [k, t) = In C(t) — n(t) Ink — 2kS(t). (11) 

Energy spectra are fitted on the intervals 2 < k < k* for 
t < 2 and on the interval 300 < k < min(/c*, fc max ) for 
t > 2, with k* = inf {k \ E{k) < 10~ 32 } denotes the be¬ 
ginning of round off noise. The fit presented in Fig.4.a in 
good agreement with the data. The time evolution of the 
fit parameter S is displayed in Fig.4.b, where it is com¬ 
pared with the exponential in time law 5(t) = 2.7e -t /°' 56 
followed by the Euler equation (see reference [24]). The 
horizontal lines show the value of the length 2/fc max 


(dashed) and 2 /k a (dotted-dashed). It is apparent that 
the model follows the Euler dynamics as long as Sk a ^ 1. 
The measure of the fit parameters is reliable as long as 
(5(f) remains larger than a few mesh sizes (<5fc max 3> 1), 
a condition required for the smallest scales to be accu¬ 
rately resolved and spectral convergence ensured. Thus 
the dimensionless quantity <5fc max is a measure of spectral 
convergence. Therefore, the self-truncation solutions are 
solution of the full partial differential equation (2) and 
not of the spectrally truncated system. 

From the mathematical point of view, Eq.(2) can be 
considered as ordinary differential equation as the right 
hand side can be shown to be a bounded continuos Lips- 
chitz map between Banach spaces. In this case, existence 
and regularity can be proved using the same techniques 
than in (finite dimensional) ordinary differential equa¬ 
tions [35]. 

D. Long-time behavior of k s t 

The self-truncation observed in Fig.2 is accompanied 
by a very slow growth of k s t until it reaches the simu¬ 
lation cut-off fc max . Such a behavior was also observed 
in the dispersive self-truncation of the truncated Gross- 
Pitaevskii equation in two [10] and three [8, 9] dimen¬ 
sions. In the Gross-Pitaevski case, the self-truncation 
wavenumber fc st was shown to obey a power-law scaling 

fcstW^f 77 , t»l- (12) 

Fig.5 displays the temporal evolution of fc st for different 
values of /3 and k a . The wavenumber k st is determined by 
the weighted average k st = \/(5/3) E a {k, t)k 2 /E a . 
With this definition, in the case of the absolute equilib¬ 
rium (8) we obtain k st = fc max . Bear in mind that k s t is 
defined for all times, even before the self-truncation starts 
to take place. It is in fact proportional to the inverse of 
the Taylor microscale [34]. 

Fig. 5.a compares k s t(t) for different values of the pe¬ 
nalization exponent /? with k a = 20 (lower three curves) 
and k a = 80 (upper three curves). For short times, when 
k s t < k a , the temporal evolution of k s t is, as expected, 
independent of (3. Furthermore, for the very short times, 
the penalization of small scales introduced by the a-term 
in the Euler-Voigt-a model is negligible and thus fc st is 
also independent of k a as the dynamics is given by the 
(standard) Euler equation (black curve with crosses). 

A behavior compatible with a power law is observed 
at long times. It is apparent that the exponent seems to 
depend on the value of /3 thought the data do not allow 
a clear determination of the exponent 77 . The power-law 
behavior is contaminated by the initial dynamic as it is 
actually expected to have the form k st (t ) ~ (t — to)' n , 
where to is the time when self-truncation starts. Simula¬ 
tions where t <0 and k s t -C fc max are difficult to reach 
with the present choice of parameters and resolution. 

In order to explore such power-law behavior a series of 
runs (19 — 21, see Table I) have been performed in resolu- 
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FIG. 5. (Color online) Time evolution of k s t(t) at resolution 
1024 3 (see Table I) . a) Different values of P for k a = 20. The 
crosses show the Euler evolution of k st (a = 0). b) p = 4 and 
different values of k a . The horizontal dashed lines represent 
the values of fc a /fe ma x. 

tion 512 3 and at small value of k a = 4. With this choice 
of parameter, self-truncation starts at low wavenumbers 
(~ k a ) and thus no Kolmogorov scaling can be observed. 
However such a choice allows very long temporal integra¬ 
tions (up to t = 2300) and to clearly observe the power- 
law behavior of k s t, as apparent in the inset of Fig.5. The 
value of the exponent clearly depends on the penalization 
exponent /?. The values measured for ry are displayed in 
Table II. The power-law observed for k s t strongly sug- 



FIG. 6 . (Color online) a) Temporal evolution of the self¬ 
similar function 'F(fe/fe s t(t)) = E a (k,t)k s t(t) (see Eq.13) for 
P = 4 and k a = 4. Data from direct numerical resolution of 
(2) at resolution 512 3 . The inset shows the temporal evolution 
of fe s t(t)/fcmax for different values of p. 

function 'F(z) is expected to behave as 'F(z) ~ z 2 for 
z <Cl and exponentially decay for z 1. Fig. 6 
displays k s t{t)E a (k,t) as a function of k/k s t(t) where 
the tendency to converge towards a self-similar distri¬ 
bution is confirmed for long times. The assumption of 
self-similarity implicitly supposes that the scale a is very 
small. Such hypothesis is valid only for long times such 
that k a -C fc s t(t) ■C fcmax- Discrepancies with the full 
self-similar form are certainly due to the finite values of 
the infra-red and ultra-violet cut-off of the simulation. 
In order to develop further this idea, in the next sections 
we introduce two theoretical models that allow to obtain 
both a clear numerical support of self-similarity and an 
analytic expression for the self-truncation exponent rj. 



cq 

II 

II 

a 

II 

Oi 

aV-Euler 

0.5 ±6 x 10 " 3 

0.25 ±3 x 10 “ 3 

0.07 ± 5 x 10 -3 

aV-EDQNM 

0.33 ±5 x 10 " 3 

0.11 ±9 x 10“ b 

0.085 ± 1 x 10 -4 

aV-Leith 

0.33 ± 10~ b 

0.15 ±2 x 10 " 4 

0.09 ±9 x 10~ b 


III. THEORETICAL MODELS 

Eddy-damped Quasi-Markovian Euler-Voigt-a 
model 


TABLE II. Values of the exponent 77 of the self-truncation 
wavenumber k s t(t) ~ t v (see ( 12 )) obtained from direct 
numerical simulation of the Euler-Voigt-a model ( 2 ), aV- 
EDQNM (14-17) and aV-Leith model (r = 2).(24). 

gests to look for self-similar behavior of the energy spec¬ 
trum E a (k,t), where the only temporal dependence of 
the spectrum is given by fc st (f). A self-similar form of 
the energy spectrum compatible with the conservation 
law (5) is given by 

Eq(M) = M*)*(m£)) (13) 

Where Eq is a constant with dimension of energy. The 


A popular model of turbulence is the so called Eddy- 
Damped Quasi-Markovian (EDQNM) closure [4, 18]. It 
was derived in the 60-70’s and is based on statistical 
closure of the velocity correlations in Fourier space plus 
some ad-hoc modeling of the dissipative time scales. It 
furnishes an integro-differential equation for the spec¬ 
trum E(k, t). EDQNM has been proved to be a powerful 
theoretical and numerical tool in the last 30 years as it al¬ 
lows to achieve a very large scale separation. It was also 
shown by Bos and Bertoglio [31] that EDQNM repro¬ 
duces well the dynamics of the truncated Euler equation, 
including the k ~ 5 / 3 and k 2 scalings together with the re¬ 
laxation to equilibrium. It was also used in [36] to give 
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an analytic prediction of the effective viscosity acting on 
the large scales of truncated Euler flows. 

The extension of EDQNM model to the Euler-Voigt-a 
case is straightforward. First, following Orszag deriva¬ 
tion [4] but using equations (2), we directly find 


dE(k, t) 
dt 


1 

1 + aPkP 


T NL (k, t) 


(14) 


where the nonlinear transfer T^l is modeled as 


T NL (k, t ) 


j j © k Pq {xy + z 3 ) 

A 


k 2 pE(p,t)E(q,t) 
1 + a^kP 


P 3 E(q, t)E{k 1 1) 
1 + aPpP 


dpdq 

pq 


(15) 


In (15), A represents a strip in (p, q) space such that the 
three wavevectors k, p, q form a triangle, x, y , z, are 
the cosine of the angles opposite to k, p, q. Q kpq is a 
characteristic time defined as 


®kpq — 


1 - exp (-fa + rj p + 7i q )t) 

Vk+Vp + Vq 


(16) 


The standard EDQNM equations are recovered by setting 
a = 0. Absolute equilibrium ( 8 ) is a stationary solution 
of (14) that satisfies detailed balance (no flux solution). 
The eddy damped inverse time rj k is defined as 


k 



s 2 E(s, t ) 
1 + a.P sP 


ds. 


(17) 


In the limit a —> 0 the standard eddy damped inverse 
time is recovered [37]. Note that for /3 —> 00 , the standard 
eddy time is also recovered for k < k a = 1/a whereas 
r)k = 0 for k > k a , consistently with the dynamics of the 
truncated Euler equation, as for k > fc max the dynamics is 
frozen. The constant A defines a time scale and we use the 
standard value A = 0.36. The truncation is imposed by 
omitting all interactions involving waves numbers larger 
than yfcmax in (15). The main difference between Eqs.(14- 
17) and the standard EDQNM model comes from the 
penalization terms that inhibit the energy transfer and 
modifies the straining. Indeed, Eqs. (14-15) can be 
rewritten as the standard EDQNM equations for the 
generalized spectrum E a by redefining the characteristic 
time (16) as Q kpq = 0 fc pg /(l+a / 3 fc / 3 )(l+a / 3 p / 3 )(l+a / 3 q' /3 ). 
The evolution of E a is thus like the one of E in the stan¬ 
dard EDQNM model but with a local modification of the 
energy transfer term. We refer to (14-17) as aV-EDQNM 
model. 

A number of simulations of the aV-EDQNM equation 
has been performed and give a behavior that is similar to 
that of the DNS of the full Euler-Voigt-a model, includ¬ 
ing comparable time scales. Fig. 7.a displays the energy 
spectrum for f3 = 2,4,6 and k a = 10 4 for a simulation 
with fc max = 43918 (and 14.2 points per octave). 




FIG. 7. (Color online) a) Energy spectra E a (k) at t = 3.5 
of the qV-EDQNM model, for different values of /3, obtained 
with k a = 10 4 and k max = 43918 (corresponding to resolution 
131072, with 14.2 points per octave), b) Temporal evolution 
of the self-similar function 'F (k/k s t(t)) = E a (k,t)k s t(t) for 
j3 = 4 and k a = 5 and kmax = 692 (28.4 points par octave). 
The inset b.l shows the temporal evolution of k a t(t) for dif¬ 
ferent values of f5. The inset b.2 shows the aV-EDQNM (red 
dashed line) and aV-Leith (solid green lines) theoretical pre¬ 
dictions (19) and (26) for the self-truncation exponent y and 
numerical data from table II. 


Fig. 7.a shows that the three zones observed in the 
Euler-Voigt-a model (see Fig.2), are also apparent in oV- 
EDQNM model spectra, with scaling laws extending for 
more than two decades. In the same spirit than in the 
Euler-Voigt-a DNS runs with a smaller value of k a = 5 
have been performed, allowing a clear determination of 
the self truncation exponent 77 . The values are presented 
in Table II. Analogously to the Euler-Voigt-a DNS, we 
look for a self-similar behavior of the energy spectrum. 
The collapse is manifest in Fig. 7.b where the self-similar 
form (13) is displayed for (3 = 4. The aV-EDQNM model 
allows for directly looking for such self similar behavior 
and find the exponent 77 by counting powers. Indeed, 
introducing the self-similar form and the variable z = 
k/k s t(t ) in (13) in the eddy damped inverse time (17) we 
































find 


riz(t) - A£’o / 2 fc st (t)^ (1 + Jk^(ty z >py dz '- (18) 

The eddy damped inverse time can thus be expressed 
as p z = XE 0 k s tl(z,ak s t)- Assuming that I(z,ak s t) ~ 
(afcst ) 7 we obtain for the characteristic time in the self¬ 
similar form Q ZkZpZq ~ '. In the same way, the non¬ 

linear transfer term (15) is found to scale with the self¬ 
truncation number as T/vx ~ k~ t 1 ~~ ,+3 ~ 2 ~ 3 ^. Equating 
the left and right hand sides of (14) we obtain k st k~ t 2 ~ 
k-'- 3i3 . Finally, the scaling k st {t) ~ t v leads to 


first contraint is given by the no-flux solution or the ab¬ 
solute equilibrium ( 8 ). Imposing that F(k) =0 for the 
absolute equilibrium E a (k ) ~ k 2 leads to m = 0 and 
n = 3/2 — q. The second contraint is for the Kolmogorov- 
like solution E a ~ fc -5 / 3 . Imposing that for such solu¬ 
tions the flux is given by F = e/(l + a^k^) r , we obtain 
p = 3/2 — q — r, being r a free parameter. Such a flux 
can be easily interpreted in the limit of large /3. Indeed, 
for k -C k a the a term is negligible, and the flux becomes 
constant. The Kolmogorov phenomenology is thus recov¬ 
ered. On the other hand for k k a the flux vanishes as 
expected. Taking into account the previous constraints, 
we obtain a family of a diffusive aV-Leith models indexed 
by the parameter r: 


1 

3/3 - 1 + 7 


(19) 


dE a 

dt 


2c/7 d 

r fc 13 / 2 d 

[ K 3/2, - 3 ] 

3 dk 

(1 + aPkP) r dk 

K 


(24) 


The exponent 7 can be computed using the ansatz 
'F(z) = z 2 exp [—z] and it reads 7 = — min [/3, 5/2]. The 
theoretical prediction confronted with the data are pre¬ 
sented in good agreement in the inset b.2 of Fig.7. 


B. a- Voigt Leith model. 


Let us now introduce another model that shares the 
same dynamics properties than Euler-Voigt-a. It is a 
spectral diffusion model that generalizes the so-called 
Leith model [19]. The original Leith model it is a phe¬ 
nomenological non-linear (local) spectral diffusion equa¬ 
tion for the energy spectrum that admits the stationary 
solutions corresponding to an absolute equilibrium and 
Kolmogorov scaling. When forcing and dissipation is 
added to the model, a steady state containing mixture 
of constant flux and thermal equilibrium was observed in 
[38]. It also known to posses self-similar solutions [39]. 

The simplest generalization of the Leith model to take 
into account the a-terrn that conserves the total energy 
E a (k ) is given by 


dE _ 1 dF 

dt (1 + a^kP) dk ’ 


( 20 ) 


where F(k ) is a spectral (non-linear) flux. Following 
Leith’s original derivation, we assume that the spectral 
flux is defined in terms of a diffusion coefficient D(k ) and 
a potential Q{k) such that 

m = (21) 

Assuming locality in the flux, by dimensional analysis 
we obtain 


D = k 9 / 2 ~ m (E / k 2 ) n (1 + a?kP) p (22) 

Q = k m {E/k 2 ) 3/2 - n {l+a f) kP) q (23) 

The dimensionless coefficient 7 sets the global time-scale 
and n, m, q,p are free parameters to be determined. The 


The standard Leith model is recovered by setting /3 = 0 
and rescaling the time. Note that Eq.(24) can be rein¬ 
terpreted as a standard nonlinear diffusion equation for 
the generalized energy spectrum E a by introducing the 
triple decay time T?(k) (see e.g. ref. [40]). In terms of 
r 3 “(fc) Eq.(24) reads 


dE a 

dt 


2 gy d 

T dk 


T«{k)k 7 E a — [k~ 2 E a 


(25) 


where the triple decay time is given by T^(k) = 

T Nh(k)/(l + aPkP) r with the eddy turnover time defined 
as r^ L (fc) = [k^^Ea 2 ]^ 1 [41]. Different choices of if 
could be used to model for instance in magnetohydrody¬ 
namics flows the presence of Alfven waves [40]. 

The solutions of the aV-Leith model (24) indeed re¬ 
produce the ones of the Euler-Voigt-a equations. Figure 
8 .a displays the energy spectrum for different values of 
/3 and r = 2. The Kolmogorov fc -5 / 3 , the equilibrium 
k 2 regimes and the fast decay for large k is manifest. 
The absence of a dissipative zone is apparent in Fig. 8 .a, 
when compared both with the aV-EDQNM case (see Fig. 
7.a) and the Euler Voigt-a model (see Fig.2). This is cer¬ 
tainly due to the locality in Fourier space of the aV-Leith 
model. 

As in the previous models, we look for self-similar so¬ 
lution of (24). Introducing (13) in to (24) we obtain in 
the limit of k k a for the self-truncation exponent 


1 = 


1 

fir — 1 


(26) 


This prediction coincides with the one of EDQNM for 
r = 2 + 7 /^ (see Eq.(19)). The self-similarity behavior 
of E a (k,t ) is apparent in Fig. 8 .b, where the self-similar 
form is displayed for /3 = 4 and r = 2. The inset shows 
the temporal evolution of k s t(t) for different values of /3. 
A power-law growth is manifest. The measured values of 
the exponent p presented inTable II are in good agree¬ 
ment with the prediction (26). 

The self-similar analysis leads to a non-linear second- 
order ordinary differential equation for 'F(z) (see (13)). 
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FIG. 8. (Color online) a) Energy spectra E a (k) of the aV- 
Leith model (24) for different values of (3 obtained with k a = 
400 and kmax = 4000. b) Temporal evolution of the self 
similar function 'F (k / k a t(t)) = E a (k,t)k s t(t) for p = 4 and 
k a = 2. The inset shows the temporal evolution of k s t(t) for 
different values of p. 


This equation cannot be solved analytically but an 
asymptotic analysis predicts ’F(z) ~ z 1 2 for 2 < 1 and 
T(z:) ~ (cte — z) 3 / 2 for z > 1. Data is compatible 
with this result (not shown). Note that, unlike the 
Euler-Voigt-a and aV-EDQNM models, the Leith model 
presents a sharp cut-off instead of an exponential decay. 


IV. CONCLUSION 

In summary the Euler-Voigt-a model allowed us to 
show that its self-truncation regime reproduces the be¬ 
havior of the truncated Euler equation [5]. We also found 
evidence for self-similarity in the long-time behavior of 


the energy spectrum. Introducing two different simpli¬ 
fied models, the aV-EDQNM model and the aV-Leith 
model, we were able to show that they present behav¬ 
iors similar to that of the Euler-Voigt-a model. We were 
able to determine the analytical values of the self-similar 
exponents of the simplified models. 

In the present work we have used only integer values 
for (3. As was noted in Sect.II A, choosing /3 = 11/3 
yields an absolute equilibrium E(Jc) ~ fc -5//3 and, in 
two-dimensions, (3 = 2/3 also yields Kolmogorov scal¬ 
ing. This can represent an interesting alternative to the 
fractal decimation method that was used in reference [22]. 

The present work can be naturally extended to the 2D 
and 3 D Ideal MHD equations. In this context, it was 
recently shown that dynamo action can be triggered by 
turbulence in absolute equilibrium [42]. Thermalization 
with (3 = 11/3 would allow for a more realistic velocity 
spectrum, mimicking the infinite Reynolds limit. 

As already stated in the introduction, sevral dissipa¬ 
tive Navier-Stokes-Voigt-a regularizations have been pro¬ 
posed as subgrid-scale models of classical turbulence clo¬ 
sure problems both in hydrodynamics [13] and in magne¬ 
tohydrodynamics [14]. We now end this paper with a few 
remarks on the relation of the present work with these 
more standard subjects. 

First, all the calculations presented above were per¬ 
formed without dissipation. This methodology clearly 
poses a challenge when one is interested in the Reynolds 
number (Re) of the flows, as e.g. in the implicit Large 
Eddie Simulations of reference [43]. However, in fully de- 
velopped turbulence [34] the dissipative scale kd is related 
to the integral scale ki by kd ~ kjRe 3 / 4 . Thus a first 
guess for a Reynolds number could be Re ~ (fc^/fc/) 4 / 3 , 
with kd given by Eq.(10) above (see section II C).This im¬ 
portant problem is beyond the present work and clearly 
deserves further attention. 

A second question is how the sweeping and straining 
processes (for review articles, see e.g. [44, 45]) need to be 
modified in the Euler-Voigt-a and related fluid models. 
This problem was discussed, for the a-Voigt EDQNM 
model, in the paragraph following Eq.(17). The identi¬ 
fication of sweeping and straining processes in the com¬ 
plete Euler-Voigt-a model Eq.(2) is an important prob¬ 
lem that must also be left for further studies. 
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